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We have studied diffusion pathways of a silyl radical adsorbed on the hydrogenated Si (100)- 
(2 x f ) surface within density-functional theory. The process is of interest for the growth of crystalline 
silicon by plasma-enhanced chemical vapor deposition (PECVD). Preliminary searches for migration 
mechanisms have been performed using metadynamics simulations. Local minima and transition 
states have been further refined by using the Nudged-Elastic-Band method. Barriers for diffusion 
from plausible adsorption sites as low as 0.35 eV have been found, but trap states have also been 
spotted, leading to a more stable configuration, with escape barriers of 0.80 eV. Diffusion through 
weakly bound physisorbed states is also possible with very low activation barriers (<50 meV). 
However, desorption mechanisms (either as SiH3 or as SIH4) from physisorbed or more strongly 
bound adsorption configurations turn out to have activation energies similar to diffusion barriers. 
Kinetic Monte Carlo simulations based on ab-initio activation energies show that the silyl radical 
diffuses at most by a few lattice spacing before desorbing at temperatures in the range 300-1000 K. 

PACS numbers: 81.15.Gh 68.47.Fg 68.43.Bc 



I. INTRODUCTION 

Plasma-enhanced chemical vapor deposition (PECVD) 
from silane is a widespread technique employed to grow 
thin films of amorphous silicon^ High growth rates are 
made possible by the deposition of reactive radicals pro- 
duced in the plasma, as opposed to conventional CVD 
where the less reactive silane is directly adsorbed at 
the growing surface. By suitably controlling the energy 
of the ions while keeping a high density of reactants, 
a related technique, called low-energy plasma-enhanced 
chemical vapor deposition (LEPECVD)^^ has intro- 
duced the possibility to obtain device-quality epitaxial 
films of crystalline silicon or silicon-germanium alloys, at 
temperatures much lower (600 °C) than those necessary 
for conventional, thermal chemical vapor deposition (800 
°C). SiH3 is supposed to be the most abundant radical 
species in the plasma discharge in most of the experimen- 
tal conditions&i. A description of the interactions of the 
silyl radical with the crystalline Si (100) surface is then of 
great interest to model the epitaxial growth at conditions 
of both PECVD and LEPECVD. 

The silyl radical can interact with the surface giving 
rise to different processes: it can remove a hydrogen 
from a saturated dimer via an Eley-Ridcal mechanism, 
or it might adsorb on the surface, where it would further 
evolve by diffusing, decomposing with hydrogen release 
at the surface or by desorbing (either as SiH3 or SfHi). 

The fate of the silyl radical depends on the degree of 
surface hydrogenation at the landing site. At low hydro- 
gen coverage the silyl can adsorb on a silicon dangling 
bond and decompose easily (into SiH2 for instance^). 
The sticking probability is supposed to decrease at higher 
hydrogen coverage. The decomposition of the adsorbed 
radical is also supposed to be strongly hindered by the 
lack of free silicon surface atoms to which H might be 
transferred. Diffusion of the adsorbed silyl from hydro- 



gen rich to hydrogen poor regions might be a conceiv- 
able route for SiH3 decomposition and insertion into the 
growing film. Fast diffusion of SiH3 adsorbed on the 
hydrogenated H:Si(100)-(2xl) surface has been recently 
predicted by first principles calculations^^. An equally 
fast diffusion of silyl radicals at the hydrogenated sur- 
face of amorphous silicon is also often invoked as the rea- 
son behind the high smoothness of the amorphous films 
grown by PECVDiS. Although several processes involv- 
ing the adsorbed silyl radical at the Si(100) surface have 
already been investigated in literature, previous theoret- 
ical studies always considered a-priory guess of reaction 
pathways, a procedure which might overlook unexpected 
mechanisms for diffusion and reaction of the adsorbed 
species^^^ii^^ 

In this paper, we investigate further the diffusion of 
SiH3 adsorbed on the H:Si(100)-(2xl) surface by making 
use of the ab-initio metadynamics tcchniqu o 17 ' 18 , a new 
simulation tool which allows for extensive search of diffu- 
sion pathways. Local minima visited during metadynam- 
ics trajectories have been then optimized and migration 
barriers between different minima further refined by the 
Nudged Elastic Band (NEB) method^ , Competitive 
mechanisms for diffusion and desorption have been inves- 
tigated by Kinetic Monte Carlo simulations which pro- 
vide the average diffusion length the silyl radical travels 
before eventually desorbing. 

After a brief description of our theoretical framework 
in Section [TTJ we present in Section IIII Al our results on 
the diffusion mechanisms for the adsorbed silyl radical, 
on the existence of trap states which hinder the radi- 
cal mobility and on desorption processes which compete 
with diffusion. In Section III we report the results of the 
Kinetic Monte Carlo simulations based on the ab-initio 
activation energies which has allowed estimating the ef- 
fect of possible errors in the calculated activated energies 
on the fate of the silyl radical. Sec. IV is devoted to our 
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FIG. 1: (color online) The slab model of the 
H : Si (100) — (2 x 1) surface, from a top a) and side b) 
views. 




conclusions. 



II. COMPUTATIONAL DETAILS 

Calculations have been performed within the frame- 
work of Density Functional Theory (DFT) with a 
gradient corrected exchange and correlation energy 
functional (PBE)2i as implemented in the codes 
PWSCF-22. for geometry optimizations and NEPJ^ cal- 
culations and CPMD22 for Car-Parrinello^ metadynam- 
ics simulation s 1 ^ 18 . Norm-conserving^ and ultrasoft2£ 
pscudopotcntials have been used for silicon and hydro- 
gen, respectively. Kohn-Sham orbitals are expanded in 
plane waves up to kinetic energy cutoff of 25 Ry. Selected 
calculations have been repeated with a norm conserving 
pseudopotential for H as well and a kinetic cutoff of 30 
Ry. 

As preliminary tests of our framework, we have cal- 
culated the reaction enthalpies for some decomposition 
reactions of disilanc (Si2Hg) within a cubic supercell of 
edge 15.3 A. Both PBE and BLYB2L22 have been used for 
these benchmark calculations on molecular systems. Ge- 
ometries have been optimized until the maximum resid- 



TABLE I: Reaction energies (AUo) and enthalpies {AH, in 
eV) calculated for some gas-phase decomposition reactions of 
disilane. Results obtained with PBE and BLYP functionals 
are compared to experimental data and previous quantum- 
chemical results. 



Si 2 H 6 


SiH 2 + SiH 4 


H2 + Si2H4 


SiH 3 + SiH 3 


AUo PBE 


2.559 


2.169 


3.147 


AUo PBE a 






4.158 


AUo BLYP 


2.213 


2.099 


2.984 


AUo MP2 6 


2.436 


2.375 


3.182 


AH PBE 


2.445 


2.033 


3.028 


AH BLYP C 


2.098 


1.963 


2.865 


AH MP2 6 


2.355 


2.064 


3.048 


AH B3LYP d 


2.259 


2.025 


3.265 


AH Exp. d 


2.355 


2.021 


3.326 



"Spin restricted calculation 
6 MP2/6-311++g** results from Ref^S 

c Vibrational contributions to enthalpy from PBE normal modes. 
d from Ref— L 



ual force on every atom was less than 0.01 eV/A. For- 
mation enthalpies have been calculated taking into ac- 
count the vibrational contribution to the quantum parti- 
tion function in the harmonic regime, with normal modes 
obtained, in turn, by diagonalization of the dynamical 
matrix, built from the numerical derivatives of the forces 
with respect to finite atomic displacements (0.015 A). 
Rotational and translational contributions to reaction en- 
thalpies have been added in the classical limit. Results 
for PBE and BLYP functionals are compared in Table U 
with experimental data and previous quantum-chemical 
calculations. The results arc in good agreement with ex- 
periments and higher levels of theory, PBE performing 
overall better than BLYP functional. For reactions in- 
volving homolithic bond breaking such as the decompo- 
sition Si2H6 — > 2SiH3, spin unrestricted calculations al- 
lowing spin polarization arc mandatory to reproduce the 
correct reaction enthalpies. Neglect of spin polarization 
(spin restricted calculations) introduces errors as large 
as 1 cV (cf. Table I) in the reaction energies. All the 
surface calculations have been then performed in a spin 
unrestricted framework (LSD-PBE). 

The H:Si(001)(2xl) surface is modeled in a slab geom- 
etry with 3D periodic boundary conditions (see FigQ]). 
The theoretical equilibrium lattice parameter obtained 
from a bulk calculation with the k-point sampling cor- 
responding to the supercell T-point is 0.7 % larger than 
the experimental one, while the theoretical value at full 
convergence in BZ integration is 0.6 % shorter than the 
experimental one. Thus, in the slab calculation we chose 
the experimental lattice parameter (5.43 A22) as an aver- 
age value of the theoretically lattice constants that can be 
obtained with different supercell sizes and k-point sam- 
pling. In a previous work^, we checked that by chang- 
ing the lattice parameter by 0.7 %, the adsorption en- 
ergy of the silyl radical on the clean Si(001)-(2xl) surface 
changes by less than 10 meV. The periodically repeated 
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FIG. 2: Minimum energy path obtained by NEB optimization 
for the migration from the minimum (a) to the minimum (b) 
of Figure [4] The energy profile with respect to a suitably de- 
fined reaction coordinate q is reported. Path A is optimized 
with the T-point only. Path B corresponds to the total ener- 
gies obtained with one special fc-point in the surface BZ for 
the geometries of path A. Path C has been fully optimized by 
using one special fc-point. 



previously visited by the trajectory in CV space, which 
acts directly on the ionic coordinates R/ as 




slabs arc separated by vacuum, 11 A wide. The slab 
has six silicon layers, each containing 12 atoms. The 
hydrogenated top surface, reconstructed in the (2x1) ge- 
ometry, contains six hydrogenated silicon dimers. The 
bottom surface of the slab is saturated by symmetric 
Siti2 groups. The S1H2 groups and the underlying silicon 
layer were kept fixed at the ideal bulk positions. Only the 
supcrccll T-point has been considered in Brillouin Zone 
sampling for geometry optimizations. We checked the 
convergence of our results with respect to Brillouin Zone 
integration by optimizing selected geometries with the 
special fc-point (in crystallographic coordinates^!) . 

Calculations with large surface supercell (three rows of 
five dimers each) have also been performed for selected 
systems. 

To uncover possible diffusion pathways of the silyl 
radical, we have made use of the metadynamics tech- 
nique which allows large barriers to be overcome in an 
affordable simulation time (few picoseconds) within ab- 
initio molecular dynamics simulations . 1 7 i 18 The method 
is based on a coarse-grained, non-Markovian dynamics in 
the manifold spanned by few reaction coordinates biased 
by a history-dependent potential which drives the system 
towards the lowest saddle point. The main assumption is 
that the reaction path could be described on the manifold 
of few collective variables (CV) 5 a ({R/}), function of the 
ionic coordinates R/ . The Lagrangian £0 of the system is 
then modified, introducing an history-dependent biasing 
potential, which affects the dynamics so as to discourage 
the system from remaining in the region already visited 
and pushes it over the lowest energy barrier towards a 
new equilibrium basin. Many variations over these basic 
principles have been explored; in this paper we use the 
straightforward approach of Refi^, which simply intro- 
duces a repulsive potential (C — Cq — g({R/},t)) built 
from the superposition of Gaussians centered at points 



g({Ri},t) = w^ex 

tj<t 

(1) 

To study the diffusion of the adsorbed silyl radical, we 
have chosen as collective variables the (x, y) surface po- 
sition of the silyl silicon atom (two CVs). The Gaussian 
hills parameters (Eq. [T]) have been chosen as As = 0.32 
a.u. and to = 0.11 cV; a new hill is added each time 
the trajectory reaches a point 2As far from the previ- 
ous Gaussian in CV space, or at worst every 250 time 
steps. We have performed Car-Parrinello metadynamics 
simulations with a time-step of 6 a.u., effective mass for 
electrons of 600 a.u., and deuterium mass for hydrogen. 
Constant temperature (300 K) on ions is enforced by a 
Nose-Hoover thermostat^ 

Although metadynamics allows computing activation 
free-energies from a finite temperature simulation, long 
simulation time (with small Gaussian height to) are 
needed to obtain accurate estimates of activation free 
energies. Actually, our metadynamics simulations are 
aimed at just obtaining a good starting guess for the 
transformation path. Then, the geometry and activation 
energy of the transition state identified along the dynam- 
ical trajectory have been further refined by using a stan- 
dard technique with fixed end points, the Nudged Elas- 
tic Band (NEB) method^ Climbing image and variable 
springs^ have been used thoroughly (but for processes 
with a barrier below 0.10 eV, where a constant spring 
constant of 0.6 a.u. has been used), with spring constants 
kmax = 0.6 a.u. and k m i n = 0.3 a.u.. A minimization 
scheme has been applied until the residual total forces 
acting on each image in direction perpendicular to the 
path was less than 0.05 cV. We have relaxed the Mini- 
mum Energy Path (MEP) using T-only calculations, then 
we have performed a self-consistent electronic structure 
optimization along the MEP using one special fc-point. 
Full optimization of the MEP geometry with one spe- 
cial point (cf. Fig. [2] for a selected process) introduces 
changes in the activation energies of the order of a few 
tens of meV with respect to the calculation with one spe- 
cial point on the gamma-point geometry 

To obtain some figures describing qualitatively the be- 
havior of the radical, and to test how much the fate of 
the radical is affected by variations of the calculated ac- 
tivation energies within the uncertainties of DFT (up 
to 0.1 eV). we have performed Kinetic Monte Carlo 
(KMC)2i simulation based on the reaction scheme from 
metadynamics and NEB calculations. 

KMC is a simulation technique which permits to at- 
tain large length and time scales, relying on the knowl- 
edge of the rates for all the relevant reaction mechanisms 
between local minima. The dynamics is performed step- 
wise, choosing at every step one of the possible mecha- 
nism with a probability proportional to its rate, and in- 
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crementing the simulation time by the appropriate time 
step following the scheme of Refill. 



III. RESULTS 
A. Mechanisms for diffusion and desorption 

The starting geometry of our mctadynamics simula- 
tions corresponds to a S1H3 radical adsorbed in the con- 
figuration of Fig. |4ja) previously identified as a pos- 
sible adsorption geometry from first principles static 
calculations^ and from direct simulations of the impact 
of the silyl on the H:Si(100)2xl surface^. In fact, we 
have shown by targeted ab-initio MD simulations that 
the S1H3 sticks on the surface in the configuration of 
Fig. HJa) once it impinges on this site with a transla- 
tional kinetic energy in the range 0.1-0.2 eV.— 

In a mctadynamics run 10 ps long, we have observed 
the diffusion of the S1H3 radical via configurations not 
previously considered in literature such as those sketched 
in Figure [3j Indeed, inspection on the molecular dynam- 
ics trajectories has allowed identifying several possible 
local minima and migration paths for the silyl radical 
which we have then refined by geometry optimization 
and NEB calculations as described below. 

Sketches of the local minima resulting from geometry 
optimization of the mctadynamics snapshots are shown 
in Fig. 2] (in the following, we will refer to "configu- 
ration (x)" as "the configuration sketched in Figure [4] 
(x)" ) along with pictures of other structures taken from 
literature, or guessed by analogy with those identified in 
the metadynamics simulation. The energy of such min- 
ima are reported in Table [Til As anticipated in section 
UTI we have checked the convergence of our results with 
respect to the size of the supercell and to the k-points 
sampling of the Brillouin Zone (BZ). Tests with more 
than one special fc-point have revealed minute changes 
in energy differences (the desorption energy from the (a) 
site changes from 0.147 eV to 0.154 eV by using one or 
four special fc-points in the total energy calculation on 
geometries optimized with T-point only). On the other 
hand, to achieve reasonable accuracy in activation ener- 
gies, it is sufficient to use one special point on geometries 
optimized with gamma-point only (cf. Figure [5]). 

The calculated activation energies for the various mi- 
gration processes among the different local minima are 
listed in Table [TTTT and a schematic (and crowded) sum- 
mary of the mechanisms examined is drawn in Figure [HI 

Large disagreement is found for the adsorption energy 
with respect to some ab-initio results previously reported 
in literature. For instance in Ref£ an adsorption energy 
of 0.63 and 0.75 eV are reported for silyl in configura- 
tions (a) and (g), while we obtain the values 0.15 and 
0.35 eV, respectively. Similarly, the energy gain for the 
removal of a hydrogen by a silyl radical is as high as 0.86 
eV in RefJ^ to be compared with our result of 0.35 eV. 
A possible source of discrepancy might be the use of a 



FIG. 3: (color online) Four snapshots from the metadynamics 
trajectory of silyl diffusion on H:Si(100)-(2xl): a) t = ps, 
b) t = 2.1 ps, c) t = 3.5 ps, d) t = 10 ps. 



a) 



b) 



c) 



d) 
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FIG. 4: (color online) Sketches of ten different local minima 
for SiH3 adsorbed on H : Si (100) — (2x1) surface. Side views 
of (a) and (i) structures are also reported at the bottom to 
illustrate the difference in geometry between strongly bound 
and physisorbed minima. The structures (e) and (g) are not 
perfectly symmetric with respect to the dimer plane, so that 
two corresponding mirror structures exist; this small asymme- 
try is due to the balance between the interaction of the SiH3 
hydrogen atoms with the surface hydrides, and the barrier 
for interconversion between the mirror configurations is likely 
to be very low. For what concerns diffusion, the two mirror 
images can be regarded as a single, symmetric structure. 




TABLE II: Energies (eV) of the local minima shown in Fig- 
ure 2] The zero of energy corresponds to SiH3 far from 
the surface, as computed from free surface and isolated silyl 
within the same supercell used to model the adsorbed silyl. 
S1H4 + db indicates the reaction energy of silane formation 
from a silyl radical abstracting a surface hydrogen atom. Col- 
umn A refers to the results that are used throughout the pa- 
per, while other columns refer to various tests we performed 
on selected structures by changing fc-points sampling and su- 
percell size. Namely, column A refers to calculations with 
6 dimers supercell and energies computed with one special 
fc-point on geometries optimized with the F point only. Col- 
umn B refers to the same calculations as in A with energies as 
well as geometries computed with the T point only. Column 
C refers to calculations with 6 dimers supercell with energies 
and geometries optimized with one special fc-point. Column 
D refers to calculations with a 15 dimers supercell, energies 
and geometries optimized with F point only. Column E refers 
to calculations as in column B but without spin polarization, 
i.e. within a spin-restricted framework (LDA-PBE instead of 
LSD-PBE). 





A 


B 


C 


D 


E 


a) 


-0.15 


-0.23 


-0.14 


-0.16 


-0.64 


b) 


-0.18 


-0.38 




-0.23 




c) 


-0.15 


-0.32 




-0.20 




d) 


-0.25 


-0.45 


-0.26 


-0.30 




e) 


-0.34 


-0.43 




-0.36 




f) 
g) 


-0.60 
-0.37 


-0.67 
-0.54 




-0.61 


-0.64 


h) 


-0.06 


-0.07 








i) 


-0.04 


-0.05 








j) 


-0.07 


-0.08 








SiH 3 


0.00 


0.00 


0.00 


0.00 


0.00 


SiH 4 + db 


-0.50 


-0.52 


-0.50 




-0.80 



spin restricted framework in Ref^^. Indeed, by repeat- 
ing our calculation with no spin polarization we have ob- 
tained adsorption and reaction energies similar to those 
reported in previous works (column E in Table [TTJ) . As 
discussed in SecHU spin unrestricted framework is, how- 
ever, recommended to properly describe an open shell 
system like the silyl radical, especially in the gas phase. 

Except for the overcoordinatcd configuration (a), the 
local minima can be classified into three categories: con- 
figurations with a broken dimer ((e), (f), (g)), and the 
unpaired electron mostly localized on a surface Si atom, 
configurations with a broken backbond ((b), (c), (d)), 
and the unpaired electron localized onto a subsuperfi- 
cial atom (see Figure [5J , and finally weakly bound "ph- 
ysisorbed" configurations of the silyl radical ((h), (i), (j)). 
Configurations in the first category appear to be lower in 
energy than those in the second one, because the sub- 
strate is less strained and the H — H repulsion is lower. 

Configurations with adsorption energies in the range 
0.15-0.35 eV are referred to hereafter as strongly bound 
adsorption states, while traps are states with still higher 
adsorption energies. The three "physisorbed" states 
could be reached from the gas phase with no activa- 
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FIG. 5: (color online) Spin polarization (m (r) = p-\ (r) — p± (r)) for the four configurations (a), (b), (f) and (j) of Figure[4] is 
shown in the four corresponding panels. Isosurfaces at m — ±1CP 3 a.u. are reported, with colors depending on the spin sign. 
In the three structures, J d 3 r |m(r)| is approximately 1.2 au, corresponding to the unpaired electron introduced by the silyl, 
plus some polarization effects. Structure (a) shows the most pronounced spin derealization, with the unpaired electron shared 
between the silyl silicon and two surface silicon atoms. 




(f) 



a) 



tion barrier. However, because of the very tiny ad- 
sorption energy, adsorption at these sites might have a 
low cross section. Indeed, these states have not been 
observed in our previous MD simulations of SiH3 im- 
pinging on the surface with energy in the range 01.-0.2 
eV— . The existence of a precursor state for silyl ad- 
sorption on hydrogenated surfaces has been proposed in 
early models of PECVD growth^, although by assum- 
ing a physisorption geometry involving a three centers 
configuration (H^Si — H — Si sur f ) which has been ruled 
out by later ab-initio calculations on H : Si(lll)ii and 
H : Si (100) 2 x 1— surfaces. The silyl can move from 
physisorbed states to configurations with a higher ad- 
sorption energy by overcoming tiny barriers (< 0.05 eV) 
which can be ascribed to hydride-hydride electrostatic 
repulsion. 

To complete the analysis of the adsorbed silyl, we have 



also considered desorption processes, either as SiH 3 or as 
SiFLj, with the concurrent removal of a surface hydro- 
gen. We have not been able to obtain a direct desorption 
mechanism from the strongly bound adsorption states, 
as paths optimized with the NEB method always show 
the physisorbed states as intermediates. Desorption as 
silyl from physisorbed configurations (h), (i) and (j) takes 
place with no barrier other than the binding energy of the 
free radical. Energy barriers of 0.07, 0.05 and 0.11 eV 
have to be overcome by the silyl radical to leave the sur- 
face as a silane molecule by removing a nearby hydrogen 
from configurations (h), (i) and (j), respectively. 

To simplify the discussion of Figure [6] which reports all 
the migration events we have investigated, at first let us 
consider jumps among strongly bound states. 

Starting from configuration (a), we can identify two 
diffusion pathways, one parallel to the dimer rows (Fig- 
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FIG. 6: (color online) A synoptic scheme of the different migration mechanisms we have considered. The energy of local minima 
and transition states (in eV) are reported with reference to the desorbed SiH3 (zero of energy). Activation energies for jumps 
in both directions can be obtained as energy differences. Transition-state energies for desorption from the physisorbed states, 
either as silyl or as silane, are also reported (large arrows). 




urc [7] (a)) and the other in the perpendicular direc- 
tion (Figure [TJd). The overall activation barrier is 
0.36 eV in both cases, corresponding to the key-event 
(b) — > (a). The direct (b)— >(b') conversion has a higher 
barrier than a mechanism which takes place via re- 
versible hydrogen exchange between the dimer atoms 
((b)— >(d)— >(d')— >(b')), which was spotted by metady- 
namics. 

However, from configuration (b), the radical can reach 
the more stable state (f) state, following one of the two 
paths highlighted in Figure El both of which have acti- 
vation energies lower than the activation barrier for dif- 
fusion. The deepest local minimum (f) acts as a "trap" 
state which hinders diffusion. The overall activation en- 
ergy for diffusion from (f) site to an equivalent (f) site 
along the dimer row is as high as 0.77 eV. Even though 
the presence of the trap state rules out fast diffusion 
at low temperatures, a barrier of 0.7 eV is still suffi- 
ciently low to be easily overcome at room or growth 
temperatures. Still, we have not yet included the ph- 



ysisorbed states into the picture. Actually, some diffu- 
sion pathways involving physisorbed states as interme- 
diates have slightly lower barriers than the mechanisms 
involving strongly bound states only, e.g. the (a) — > (a') 
path shown in Figure with an overall barrier of 0.17 
eV. 

The evolution of the adsorbed SiH3 radical we have dis- 
cussed so far can be summarized by the simplified scheme 
of Figure [10] A high barrier has to be overcome to escape 
from the realm of trap states and reach some "mobile" 
local minima. From one of those states, with similar bar- 
riers, the system might fall back into the traps, diffuse 
to nearby strongly bound states or reach one of the ph- 
ysisorbed states. Barriers for migration between different 
strongly bound states (~ 0.3 eV) are similar to barriers 
for migration between strongly bound and physisorbed 
states. From the latter physisorbed states, with barriers 
below 0.1 eV, the silyl can go back to the strongly bound 
states, diffuse amongst the physisorbed states or desorb, 
with similar barriers, cither as silyl or as silane. The 
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FIG. 7: (color online) Diffusion pathways for a SiH3 radical, a) parallel to the dimer row, and b) perpendicular to the dimer 
row. Energies of transition states and intermediate minima (in eV) are given. The zero of energy corresponds to the desorbed 
SiH 3 . 



0,17 -0.03 -0.02 -0.03 0.17 




FIG. 8: (color online) Two paths leading to the trap state 
(f) from the local minimum (b). The zero of energy (eV) 
corresponds to the desorbed S1H3. 




accuracy of our calculations is not sufficient to identify 
which are the dominant processes at low temperature, 
since uncertainties of the order of 0.1 eV in the activa- 
tion energies are expected due to approximations in the 
exchange and correlation functional. However, at typi- 



FIG. 9: (color online) A path for diffusion between two equiv- 
alent (a) sites, passing through physisorbed state (i). Energies 
(eV) of minima and transition states are given. The zero of 
energy corresponds to the desorbed S1H3. See Figure [4] for 
sketches of the two structures viewed along [Oil]. 

0,02 ODO 0,02 




cal growth temperatures of PECVD and LEPECVD the 
different events are likely to have similar probability. To 
estimate the average diffusion length that a silyl radical 
can travel before desorption takes place, we have per- 
formed kinetic Monte Carlo simulations using the acti- 
vation energies for the different processes summarized in 
Table III. The results are discussed in the next section. 

As a final remark, we compare our results on diffusion 
mechanisms with previous ab-initio works appeared in 
literature. The migration (a) — ► (a') between two adja- 
cent dimers rows and the (g) — > (g') jump within an open 
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TABLE III: Activation energies (in eV) for various migration 
mechanisms of SiH3 on H : Si (100) — (2x1). Column B re- 
ports values obtained by NEB optimization of the minimum 
energy path by sampling the BZ at the F-point only. Column 
A corresponds to the activation energies computed with one 
special k-point on geometries optimized with T-point only. 
The primed and not-primed configurations are equivalent by 
symmetry: (a) — * (a') corresponds to a jump across the chan- 
nel, (b) — ► (b') and (d) — > (d') corresponds to the conversion 
between two structures equivalent with respect to the mir- 
ror plane orthogonal to [Oil], (g) — * (g') corresponds to the 
migration of the radical between the two atoms of the bro- 
ken dimer, (i) — * (i') and (i) — » (i") correspond to two steps 
of a jump along the dimer row, between two equivalent ph- 
ysisorbed states. 
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FIG. 10: (color online) Simplified scheme that accounts for 
most of the features of the complex reactions scheme of Fig- 
ure[6] the activation energies have to be considered just quali- 
tative estimates, because of both the limited accuracy of DFT 
calculations and the fact that each of the barriers in this fig- 
ure corresponds to different pathways with slightly different 
activation energies (cf. Figure For simplicity we label as 
" chemisorbed" (here and in Table IIV|I the bound states of 
S1H3 with adsorption energy in the range 0.15-0.35 eV. 
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TABLE IV: Energy (eV) of "minima" and "transition states" 
used in the simplified Kinetic Monte Carlo simulation based 
on Figure 1101 The zero of energy corresponds to the des- 
orb silyl ("silyl" in the table). For simplicity we label as 
"chemisorbed" (here and in Fig. I10[) the bound states of SiH3 
with adsorption energy in the range 0.15-0.35 eV. Prefactors 
for the "fixed-prefactor" model have been set equal to 1 THz, 
or twice as large in the presence of two symmetry-equivalent 
paths. We have also considered rate prefactors v* obtained 
with harmonic transition state theory for selected processes, 
namely (a)— >(b) for the reactions starting in a strongly bound 
state ("chemi" or "trap" in the table), (i)— »(i') for reactions 
starting in a physisorbed state, and (i)— > SiH3 for desorption. 
The "forward" and "backward" values for v* (THz) used in 
the simulations of Figure [13] are given. 
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dimer have been investigated in Refi^ within a framework 
very similar to ours which, however, produced activation 
energies as high as 0.9 and 0.6 eV for the two process, 
respectively, as opposed to our values of 0.3 and 0.22 eV. 
The discrepancy is due to a very different geometry of 
the transition state which in Ref.— corresponds (both for 
(a) — > (a') and (g) — > (g') jumps) to a silyl nearly flat, 
sp 2 -likc followed by an "umbrella flip" of the hydrogen 
atoms of the silyl. In our case, by choosing the trajec- 



tory obtained from the metadynamics simulation as ini- 
tial MEP in the NEB optimization, we have been able 
to find a lower transition state involving small defor- 
mation of the SiH3 radical which undergoes a rotation 
around its C3 axis during the jump still keeping a sp 3 - 
like conformation. Indeed, within the NEB method, it 
is possible that an unappropriate choice of the initial 
MEP might lead to a reaction path with an activation 
energy higher than the lowest energy path to the prod- 
ucts. Actually, we remark that an even lower barrier 
(0.17 vs 0.30 eV) for the jump across the channel be- 
tween two neighboring rows (a) — > (a') can be obtained 
along the path (a) — > (i) — > (a'), with the physisorbed 
state (i) as an intermediate. Another issue worth being 
mentioned is a discrepancy we have found with previous 
workii on the diffusion mechanism of the silyl along the 
dimer row. Using LSD-PBE (spin unrestricted calcula- 
tions) we have not been able to optimize the intermediate 
state described as a local minimum in Figure 3 of Refi^ (a 
silyl adsorbed on the five-fold coordinated Si atom of the 
Si-H group) which has been obtained probably within a 
spin-restricted (LDA-PBE) framework. In fact, the con- 
figuration proposed as a local minimum along the diffu- 
sion path in Ref 1^, transforms in configuration (e) , in our 
NEB optimization. Moreover, by searching for a direct 
(a)— >(e) transformation, we have always found configu- 
ration (b) as an intermediate minimum. Consequently, 
we have not been able to find a minimum energy path 
for a "direct" jump (a) — ► (a") along the dimer row as 
reported in Refi^. Probably the discrepancies with pre- 
vious work might still be ascribed to the neglect of spin 
polarization in RefjS,. 
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B. Kinetic Monte Carlo 

Based on the reaction scheme of Figure [6] we have per- 
formed Kinetic Monte Carlo (KMC) simulations of dif- 
fusion and desorption of the adsorbed silyl. We have 
considered a lattice model, each lattice site correspond- 
ing to a 2 x 1 unit cell. At each site different configu- 
rations for the silyl, corresponding to the different min- 
ima of Table 11 (and Fig. [5]) are possible. We have 
only omitted configuration (e), as it is only a shallow 
minimum along the (b) — ► (f) migration pathway, caus- 
ing frequent oscillation of the system between minima (f ) 
and (e) which slows down dramatically simulations at low 
(room) temperatures. We have estimated reaction rates 
v at different temperatures within transition state the- 
ory, i.e. v = v* exp (— E a /kBT), with activation energies 
E a given in Table III and prefactor v* set to 1 THz for all 
processes. In some cases, for example for the (a) — > (a') 
jump, two pathways equivalent by symmetry exists: we 
have used a 2 THz prefactor in this and similar cases. The 
choice of a common prefactor is obviously questionable. 
We have calculated within the harmonic approximation 
to transition state theory the prefactors for a diffusion 
process among strongly bound states ((a)— >(b)), for the 
diffusion between two physisorbed states ((i)— >(i')) and 
for desorption ((i)— >SiHa). The prefactor v* is given 

by W v j IW v j i where v~p and vf S ^ arc the pos- 
itive phonon frequencies for the initial and transition 
states respectively, obtained by diagonalizing the dynam- 
ical matrix, calculated in turns by finite displacements of 
the atoms. For the desorption process, which does not 
display a saddle point on the potential energy surface 
along the desorption path, we have instead used the vi- 
brational frequency of the phonon whose displacement 
pattern is parallel to the minimum energy path. The re- 
sults, 28 THz ((a)->(b)), 2 THz ((i)^(i')) and 0.8 THz 
((i)^SiH3), differ less than the typical errors in this kind 
of calculations, which are at least one or two orders of 
magnitude. Therefore, given the considerable compu- 
tational resources needed to compute all prefactors, we 
have judged not worthy to pursue this task at the mo- 
ment. 

Moreover, we have considered a single silyl adsorbed on 
a perfect H:Si(100)-(2xl) surface which is obviously very 
far from the real growth conditions in PECVD where 
partial surface hydrogenation and large coverage of ad- 
sorbed species are expected. Still, these simplified KMC 
simulations can provide crucial information on SiH3 dif- 
fusion to model the behavior of the silyl radical at the 
more complex conditions of real PECVD growth. 

Within our scheme, the final fate of the adsorbed silyl 
is always desorption, either as a silyl or a silane. The 
KMC simulations let us estimate the average resident 
time of the silyl on the surface and the average maxi- 
mum distance the silyl can travel before desorbing. We 
have checked the stability of the results versus random 
variations of the barriers within the typical DFT uncer- 
tainties. To this aim, we have performed several sets of 



simulations adding a uniform noise with zero mean and 
±0.1 eV width to all the minima and transition states 
energies in Figure [51 discarding the set of randomized ac- 
tivation energies whenever one of them is negative. Some 
other sets have also been discarded since they lead the 
system to oscillate for a long time (> 10 6 steps) between 
two minima separated by tiny barriers. In all simula- 
tions the silyl starts from site (a) since it corresponds to 
the adsorption site observed in dynamical Car-Parrinello 
simulations of SiH3 impinging on the surface (cf. Section 
IIIA). 

The calculated distribution of resident times (times be- 
fore desorption) is reported in Fig. QTJ as obtained by 
averaging over 10 6 KMC simulations. Results with ac- 
tivation energies from Table III and with the energies 
randomized in every simulation are compared. The dis- 
tribution of resident times, h(x), is reported as a function 
of the logarithm of the resident time t (x = logiot), i.e. 
J x ^ h (x) dx returns the fraction of trajectories with a res- 
ident time between 10 Xl and lO^ 2 . In this representation, 
for an exponential decay ke~ kt , the histogram on a log- 
arithmic ir-axis has a peak at — log 10 k. Different peaks 
on the histogram are therefore related to processes with 
different desorption rates. The distribution of resident 
times reveals three peaks at 500 K: a nearly instanta- 
neous desorption from site (a) (low t peak), a desorption 
from site (g) easily accessible from (a) (peak at interme- 
diate t) and finally desorption from trap states at longer 
t. The peak due to desorption from site (g) reduces to a 
shoulder of the first peak in simulations with randomized 
activation energies. At higher temperatures the proba- 
bility to reach the traps increases and the importance of 
longer resident times (area of the peak at the longest t) 
increases. Therefore, most of the radicals desorb quickly, 
without reaching the trap states, but a fraction of the 
adsorbed SiH3 get to the trap, and remain adsorbed for 
longer times, of the order of 10 -5 s at 500 K. The abso- 
lute position of the peaks in h{x) depends on the choice 
of the prefactor v* in the reaction rates, but the ratio of 
different positions and the relative area of the peaks do 
not. 

In Figure [T2"l we report the results on the average max- 
imum distance the radical can travel before desorbing. 
For simplicity, only displacement along [Oil] are consid- 
ered. Diffusion distances have been calculated as the 
displacement along [011] from the initial site (a); each 
minimum within the unit cell corresponds to a displace- 
ment equal to the relative coordinate of the Si atom with 
respect to site (a). Averaging over random changes of 
the barriers results in a variation of the relative impor- 
tance of diffusion and desorption processes which is sig- 
nificative at room temperature, but negligible at higher 
temperatures. This is easy to understand, as, at room 
temperature, the radical would travel very long distances 
if the barrier for desorption was even a few tens of meV 
higher than that for diffusion among physisorbed states. 
Anyway, it is very unlikely that the radical could jump 
further than few cells before desorbing. The average dif- 
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FIG. 11: Distribution of resident times t of the silyl on the sur- 
face (before desorption) . The initial configuration is always 
site (a) (see text). The results are averaged over 10 6 indepen- 
dent KMC runs. The distribution is shown with a logarithmic 
scale on abscissa {h(x),x — logiot). Panel a) and b) corre- 
spond to simulations at T = 500 K and T — 1000 K, respec- 
tively. Curve A corresponds to the full reaction scheme, with 
the energies given in Figure [6] Curve B corresponds to the 
complete scheme, where for each simulation energies for the 
minima and the transition states have been randomized (see 
text). Curve C and D corresponds to the simplified model, 
with parameters given in Table IIVI and or with randomized 
energies (see text), respectively. 
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fusion lengths in Fig. [12] are independent from the value 
chosen for the common rate prefactor v* . 

To investigate the dependence of the diffusion length 
on the variation of rate prefactors we have considered a 
simplified model corresponding to Figure [TO] with the set 
of parameters reported in Table [TV] Although the model 
is much simpler than the complete reaction scheme, there 
are still too many parameters to attempt a significative 
fitting procedure. We have thus chosen the values in Ta- 
ble IIVI hcuristically, as we only wish to show that the 
complete scheme is very redundant, and the simplified 
scheme provides a reasonable description of the behavior 
of the silyl. In the simplified model only jumps between 
adjacent cells are considered. Displacement of the silyl 
among different minima within the same unit cell are not 
included in the calculation of the maximum average dis- 
placement. Firstly, we set all the prefactors equal, as 
done in the fully detailed KMC, for sake of comparison. 
The results are reported in Figs. [TT] and [T3J The sim- 



FIG. 12: Maximum distance along [011] reached by the silyl 
before desorption as a function of temperature. The results 
are averaged over 10 6 independent KMC runs. Curve A cor- 
responds to the full reaction scheme, with the energies given 
in Figure [S] Curve B corresponds to the complete scheme 
with randomized activation energies (see text). Curve C and 
D corresponds to the simplified model, with parameters given 
in Table ITVl and or with randomized energies (see text), re- 
spectively. Distances are given in units of the lattice spacing 
along [011] (a y =3.84 A). 
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FIG. 13: Comparison of KMC simulations performed on the 
simplified model and constant prefactors (continuous curve 
A) and with prefactors chosen according to harmonic transi- 
tion state theory predictions (dotted curve B, cf. Table ITV)l . 
In both cases averages are performed over 10 6 independent 
KMC runs, with activation barriers randomized as described 
in the text, a) Distribution of resident times of the silyl on the 
surface at 500 K (cf. Figure [TT]). b) Maximum distance along 
[011] reached by the silyl before desorption as a function of 
temperature (cf. Figure I12[) . Distances are given in units of 
the lattice spacing along [011] (a y =3.84 A). 
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plificd model is able to reproduce the main features of 
the complete KMC results. Both the shape of h(x) (but 
for the smaller structure due to desorption from (g)) and 
the order of magnitude of the average diffusion length 
are well reproduced. The simplified model can thus be 
considered as a good starting point to develop a more 
elaborate KMC model for realistic conditions with vari- 
able H and SiH3 coverage. Then, we have investigated 
the dependence of the results on the choice of different 
prefactors for different reactions. To this aim we have 
used the prefactors obtained by harmonic transition state 
theory for the three representative reactions ( Table [TV)) . 
The residence-time histogram (Figure [13]) is only slightly 
shifted. The attempt frequency for diffusion is in this 
case almost three times the one for desorption, but the 
average diffusion length is only doubled. This suggests 
that the model is also relatively stable against changes in 
the estimated prefactors, making our predictions reliable, 
despite the unavoidable errors on energies and overall 
reaction rates. We can therefore conclude that desorp- 
tion events prevents diffusion of SiH 3 radicals for more 
than a few lattice spacings. Silyl migration can account 
for local rearrangement effects, resulting in an increase 
in the capture area of surface dangling bonds, but can- 
not be invoked to explain long-range diffusion or surface 
smoothening on a mesoscopic scale. 

IV. CONCLUSIONS 

We have investigated the diffusion and desorption 
mechanisms of the S1H3 radical adsorbed on H:Si(100)- 
(2x1) by means of ab-initio calculations. Preliminary 
metadynamics Car-Parrincllo simulations aided us in 
identifying local minima and diffusion pathways of the 
the adsorbed silyl. Activation energies and minimum cn- 
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